RESULTS NOTIME excluding/regrouping rare species

Running models with different datasets:

Using only the 5x5 m quadrat scale

#Figures colors/labels
cores <- c("#c00000", "#b666d2", "#70ad47", "#ffc000")

labvit <- as_labeller(c(grow = "Growth", mort="Mortality", 
                        rec="Recruitment"))
grplab = as_labeller(c(quadrat = "Space", `quadrat:sp`= "Space X Species",
                       sp="Species", Residual = "residual"))
grpvit <- as_labeller(c(grow = "Growth", mort="Mortality", 
                        rec="Recruitment",
                        quadrat = "Space", `quadrat:sp`= "Species X Space",
                        sp="Species", Residual = "Residual"))

Data

Growth

grp <- c("all","exclude", "regroup")
tudo <- list()
path1 = here("models_outputs/models_notime/")
for (i in 1:length(grp)){
  path <- paste0(path1, grp[i], "/grow/table")
  files = list.files(path)
  all <-  lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
  names(all) <- substr(files,1,nchar(files)-6)
  tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")

#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$fplot,tes$data)/5
##       
##        all exclude regroup
##   ama    1       1       1
##   bci    6       6       6
##   edo    2       2       2
##   fus    3       3       3
##   idc    1       1       1
##   kor    1       1       1
##   lam    3       3       3
##   ldw    1       1       1
##   len    2       2       2
##   lpl    1       1       1
##   luq    3       3       3
##   mos    2       2       2
##   pas    4       4       4
##   scbi   1       1       1
##   serc   1       1       1
##   sin    2       2       2
##   ucsc   2       2       2
##   wab    2       2       2
##   wfdp   1       1       1
##   wyw    3       3       3
##   zof    1       1       1
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"

grow <- tes %>% filter(term != "intercept") %>%
  left_join(intercepts[,c(1:4,9)], by= c("id", "data", "fplot", "time")) %>%
  group_by(id, data, fplot, time) %>% 
  mutate(VPC = variance/sum(variance))  %>% ungroup() %>%
  mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) %>%
  as.data.frame()

# divide the sd of the terms by the mean growth (intercept) for each dataset.
grow$stdev = grow$Estimate/grow$intercept

Mort

grp <- c("all","exclude", "regroup")
tudo <- list()
path1 = here("models_outputs/models_notime/")
for (i in 1:length(grp)){
  path <- paste0(path1,grp[i], "/mort/table")
  files = list.files(path)
  all <-  lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
  names(all) <- substr(files,1,nchar(files)-6)
  tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")

#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$fplot,tes$data)/5
##       
##        all exclude regroup
##   ama    1       1       1
##   bci    7       7       7
##   edo    2       2       2
##   fus    3       3       3
##   idc    1       1       1
##   kor    1       1       1
##   lam    3       3       3
##   ldw    1       1       1
##   len    2       2       2
##   lpl    1       1       1
##   luq    3       3       3
##   mos    2       2       2
##   pas    5       5       5
##   scbi   1       1       1
##   serc   1       1       1
##   sin    2       2       2
##   ucsc   2       2       2
##   wab    2       2       2
##   wfdp   1       1       1
##   wyw    3       3       3
##   zof    1       1       1
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"

mort <- tes %>% filter(term != "intercept") %>%
  left_join(intercepts[,c(1:4,9)], by= c("id", "data", "fplot", "time")) %>%
  group_by(id, data, fplot, time) %>% 
  mutate(VPC = variance/sum(variance))  %>% ungroup() %>%
  mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) %>%
  as.data.frame()
mort$stdev = mort$Estimate

Rec

grp <- c("all","exclude", "regroup")
tudo <- list()
path1 = here("models_outputs/models_notime/")
for (i in 1:length(grp)){
  path <- paste0(path1,grp[i], "/rec/table")
  files = list.files(path)
  all <-  lapply(files, function(x) get(assign(x, load(paste0(path, "/", x)))))
  names(all) <- substr(files,1,nchar(files)-6)
  tudo[[i]] <- all
}
names(tudo) <- grp
tes <- map(tudo,bind_rows, .id="id") %>% bind_rows() %>% filter(q_size=="quad_5")

#tes %>% count(data,fplot) %>% mutate(n=n/8)
table(tes$fplot,tes$data)/5
##       
##        all exclude regroup
##   bci    7       7       7
##   edo    2       2       2
##   fus    3       3       3
##   idc    1       1       1
##   kor    1       1       1
##   lam    3       3       3
##   ldw    1       1       1
##   len    2       2       2
##   lpl    1       1       1
##   luq    3       3       3
##   mos    2       2       2
##   pas    5       5       5
##   scbi   1       1       1
##   serc   1       1       1
##   sin    2       2       2
##   ucsc   2       2       2
##   wab    2       2       2
##   wfdp   1       1       1
##   wyw    3       3       3
##   zof    1       1       1
#remove intercept to include again as a separate column
intercepts <- tes %>% filter(term == "intercept")
names(intercepts)[9] = "intercept"

rec <- tes %>% filter(term != "intercept") %>%
  left_join(intercepts[,c(1:4,9)], by= c("id", "data", "fplot", "time")) %>%
  group_by(id, data, fplot, time) %>% 
  mutate(VPC = variance/sum(variance))  %>% ungroup() %>%
  mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual")) %>%
  as.data.frame()
rec$stdev = rec$Estimate

Excluding recruitment for Wytham Woods, as in “all data”.

rec <- rec %>% filter(fplot!="wyw")

Combining results.

comb <- bind_rows(list(grow=grow, mort=mort, rec=rec), .id="vital") %>%
  mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp","Residual")) %>%
  ungroup()

Comparing groups

comb %>% 
  ggplot(aes(x=data, y=VPC)) +
  geom_col(aes(fill=term, color=term), alpha=0.7) +
  scale_fill_manual(values=cores) +
  scale_color_manual(values=cores) +
  facet_grid(vital~fplot+time) +
  theme(axis.text.x = element_text(angle=45, hjust=1))

mean over fplots

comb %>% group_by(vital,data,term) %>% summarise(VPC=mean(VPC)) %>%
  ggplot(aes(x=data, y=VPC)) +geom_col(aes(fill=term, color=term), alpha=0.7) +
  scale_fill_manual(values=cores) + 
  scale_color_manual(values=cores) +
  facet_grid(~vital)+
  theme(axis.text.x = element_text(angle=45, hjust=1))+
  labs(tag="a)")

#ggsave("figs/FIG_S4.4a_groups_mean.jpg", width=20, height = 10, units = "cm")
comb %>%
  ggplot(aes(x=data,y=VPC, col=data)) + geom_boxplot()+
  facet_grid(vital~term) +
  theme(axis.text.x = element_text(angle=45,hjust=1),
        legend.position = "none")

comb %>% mutate(xis = as.numeric(as.factor(data))) %>%
  ggplot(aes(x=xis,y=stdev, col=paste(fplot,time))) + geom_line() +
  facet_grid(vital~term) + 
  scale_x_continuous(breaks = 1:3, labels=levels(comb$data),
                     name="") +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none")

comb %>% mutate(xis = as.numeric(as.factor(data))) %>%
  ggplot(aes(x=xis,y=VPC, col=paste(fplot,time))) + geom_line() +
  facet_grid(vital~term) + 
  scale_x_continuous(breaks = 1:3, labels=levels(comb$data),
                     name="") +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        legend.position = "none")

Difference to all data

Absolute difference

dif <- comb %>% 
  pivot_wider(id_cols = c(vital,fplot,time,term),
              names_from = data,
              values_from = c(stdev,VPC, richness)) %>%
  mutate(stdev.dif_exclude = stdev_exclude - stdev_all,
         VPC.dif_exclude =  VPC_exclude - VPC_all,
         stdev.dif_regroup = stdev_regroup - stdev_all,
         VPC.dif_regroup =  VPC_regroup - VPC_all) %>%
  select(vital,fplot,time,term,richness_exclude,
         stdev.dif_exclude,
         stdev.dif_regroup, VPC.dif_exclude, VPC.dif_regroup) %>%
  pivot_longer(6:9, names_to = c("index", "data"), names_sep = "_",
                #names_pattern = "(.*_).(.*)",
                values_to =  "value") %>%
  pivot_wider(names_from = index, values_from = value)
dif %>%
  ggplot(aes(x=data, y=stdev.dif)) + 
  facet_grid(vital~term, scales= "free_x") +
  geom_boxplot() +
  geom_quasirandom(aes(col=fplot)) +
  ylab("SD no.rare - SD all") + xlab("") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Difference in SD") +
  theme(legend.position = "right",
        panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45, hjust=1))

dif %>%
  ggplot(aes(x=data, y=VPC.dif)) + 
  facet_grid(vital~term, scales= "free_x") +
  geom_boxplot() +
  geom_quasirandom(aes(col=fplot)) +
  ylab("VPC no.rare - VPC all") + xlab("") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Difference in VPC") +
  theme(legend.position = "right",
        panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45, hjust=1))

Relative difference SD

difr <- comb %>% 
  pivot_wider(id_cols = c(vital,fplot,time,term),
              names_from = data,
              values_from = c(stdev,VPC, richness)) %>%
  mutate(stdev.dif_exclude = (stdev_exclude - stdev_all)*100/stdev_all,
         stdev.dif_regroup = (stdev_regroup - stdev_all)*100/stdev_all) %>%
  select(vital,fplot,time,term,richness_exclude,
         stdev.dif_exclude,stdev.dif_regroup) %>%
  pivot_longer(6:7, names_to = "data",
                values_to =  "stdev.dif") %>%
  mutate(data = substr(data,11, nchar(data)))
difr %>%
  ggplot(aes(x=data, y=stdev.dif)) + 
  facet_grid(vital~term, scales= "free_x") +
  geom_boxplot() +
  geom_quasirandom(aes(col=fplot)) +
  ylab("Proportional difference (%)") + xlab("") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Relative difference in SD") +
  theme(legend.position = "right",
        panel.background = element_rect(colour="lightgray"),
        axis.text.x = element_text(angle=45, hjust=1))

difr %>% group_by(vital, term, data) %>% summarise(mean=mean(stdev.dif, na.rm=T)) %>%
  pivot_wider(names_from = term, values_from = mean)
## # A tibble: 6 × 6
## # Groups:   vital [3]
##   vital data    quadrat `quadrat:sp`     sp Residual
##   <chr> <chr>     <dbl>        <dbl>  <dbl>    <dbl>
## 1 grow  exclude 11.0           8.47  -11.8      13.2
## 2 grow  regroup  8.71          3.97  -12.7      13.9
## 3 mort  exclude  0.0278       -1.03  -14.9       0  
## 4 mort  regroup -2.62         -0.591 -16.4       0  
## 5 rec   exclude -0.565        -0.299  -9.26      0  
## 6 rec   regroup  0.512        -0.768 -10.4       0

table

Relative differences in SD

difr %>% group_by(vital, term, data) %>% 
  summarise(mean=mean(stdev.dif, na.rm=T)) %>%
  pivot_wider(names_from = c(vital,data), values_from = mean) %>%
  kable(digits=1)
term grow_exclude grow_regroup mort_exclude mort_regroup rec_exclude rec_regroup
quadrat 11.0 8.7 0.0 -2.6 -0.6 0.5
quadrat:sp 8.5 4.0 -1.0 -0.6 -0.3 -0.8
sp -11.8 -12.7 -14.9 -16.4 -9.3 -10.4
Residual 13.2 13.9 0.0 0.0 0.0 0.0

Mean absolute differences in VPC

dif %>% group_by(vital, term, data) %>% 
  summarise(mean=mean(VPC.dif, na.rm=T)) %>%
  pivot_wider(names_from = c(vital,data), values_from = mean) -> write
kable(write)
term grow_exclude grow_regroup mort_exclude mort_regroup rec_exclude rec_regroup
quadrat 0.0048918 0.0039062 0.0081885 0.0061325 0.0069821 0.0130636
quadrat:sp 0.0079957 0.0028262 0.0139204 0.0162427 0.0076814 0.0044177
sp -0.0831334 -0.0908756 -0.0704275 -0.0754252 -0.0321646 -0.0368164
Residual 0.0702459 0.0841432 0.0483186 0.0530499 0.0175012 0.0193351

Comparing with NUMBER of rare species

dif %>%
  ggplot(aes(x=richness_exclude, y=stdev.dif, shape=data)) + 
  geom_point(aes(col=fplot)) + 
  geom_smooth(method="lm", se=F, aes(linetype=data)) +
  facet_grid(vital~term) +
  ylab("Difference in SD to original data") +
  xlab("log N of rare species") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Difference in SD") +
  scale_x_log10()+
  theme(axis.text.x = element_text(angle=45, hjust=1),
        )

difr %>%
  ggplot(aes(x=richness_exclude, y=stdev.dif, shape=data)) + 
  geom_point(aes(col=fplot)) + 
  geom_smooth(method="lm", se=F, aes(linetype=data)) +
  facet_grid(vital~term) +
  ylab("Relative difference in SD to original data") +
  xlab("log N of rare species") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  ggtitle("Proportional difference in SD to original data") +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        panel.background = element_rect(fill="white", color="black"))

pm <- dif %>%
  ggplot(aes(x=richness_exclude, y=VPC.dif, shape=data)) + 
  geom_point(aes(col=fplot)) + 
  facet_grid(vital~term) +
   geom_smooth(method="lm", se=F, aes(linetype=data)) +
  ylab(" Absolute difference in VPC to original data") +
   xlab("log N of rare species") +
  geom_hline(yintercept = 0, linetype= "dotted")  +
  theme(axis.text.x = element_text(angle=45, hjust=1),
        panel.background = element_rect(fill="white", color="black"))
pm

Dirichlet regression for exclude results with species richness

Richness by rarefaction at the 6ha plot size (sampling unit 20x20m)

load(here("data", "rarefaction.curves.Rdata"))
rich.rare <- rare[rare$sites ==150,2:3]
colnames(rich.rare)[1] <- "richness.rarefaction"
comb <- comb %>% left_join(rich.rare, by="fplot")

Latitude

load(here("data", "plots_structure.Rdata"))
comb <- comb %>% left_join(plots.structure[,1:2], "fplot")
#mean for each forest
mcomb <- comb %>% group_by(vital,data, fplot,term, lat) %>%
  summarise(VPC = mean(VPC),
            stdev = mean(stdev),
            richness.rarefaction = mean(richness.rarefaction))

GROWTH

predirig <-  mcomb %>% 
  filter(vital == "grow", data == "exclude") %>%
  select(fplot, term, VPC, richness.rarefaction,lat) 
diridatag <- predirig %>%
  pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
  mutate(log.rich = log(richness.rarefaction),
         log.rich.o=log.rich) %>%
  mutate_at(vars( log.rich), scale)

vpcg <- DR_data(diridatag[,c("quadrat", "sp", "quadrat:sp",  "Residual")])
#plot(vpc)
grow2 <- DirichReg(vpcg~log.rich, data=diridatag, model="alternative", base=4)
summary(grow2)
## Call:
## DirichReg(formula = vpcg ~ log.rich, data = diridatag, model = "alternative",
## base = 4)
## 
## Standardized Residuals:
##                 Min       1Q   Median       3Q     Max
## quadrat     -0.9058  -0.4992  -0.2701  -0.0166  0.9614
## sp          -1.5134  -0.7438  -0.3542   0.2168  3.7107
## quadrat:sp  -1.7533  -0.5644  -0.2376   0.3694  2.6473
## Residual    -2.0096  -0.8659   0.3287   1.1003  1.8630
## 
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.6205     0.1744 -15.026  < 2e-16 ***
## log.rich      0.4740     0.1675   2.829  0.00466 ** 
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.09746    0.09944 -11.036   <2e-16 ***
## log.rich    -0.20184    0.10053  -2.008   0.0447 *  
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.44533    0.11231 -12.869   <2e-16 ***
## log.rich    -0.08056    0.11309  -0.712    0.476    
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
## 
## PRECISION MODEL:
## ------------------------------------------------------------------
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.4105     0.1789   19.06   <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-likelihood: 102.9 on 7 df (38 BFGS + 1 NR Iterations)
## AIC: -191.7, BIC: -184.4
## Number of Observations: 21
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatag$log.rich),
                                      max(diridatag$log.rich), length.out=10))
pred <- predict(grow2, newdata = newdata, se=T)
confint(grow2)
## 
## 95% Confidence Intervals (original form)
## 
## - Beta-Parameters:
## Variable: quadrat
##                2.5%    Est.   97.5%
## (Intercept)  -2.962  -2.621  -2.279
## log.rich      0.146   0.474   0.802
## 
## Variable: sp
##                2.5%    Est.   97.5%
## (Intercept)  -1.292  -1.097  -0.903
## log.rich     -0.399  -0.202  -0.005
## 
## Variable: quadrat:sp
##                2.5%    Est.   97.5%
## (Intercept)  -1.665  -1.445  -1.225
## log.rich     -0.302  -0.081   0.141
## 
## Variable: Residual
##   variable omitted
## 
## - Gamma-Parameters
##              2.5%  Est.  97.5%
## (Intercept)  3.06  3.41   3.76
colnames(pred) <- colnames(vpcg)
newdata$log.rich.o <- newdata$log.rich*sd(diridatag$log.rich.o) +
  mean(diridatag$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredg <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")

newpredg %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
                                       "Residual")) %>%
  ggplot(aes(x=rich.o, y=pred, col=term)) +
  geom_line() + facet_grid(~term) +
  geom_point(data=predirig, aes(x=richness.rarefaction, y=VPC, col=term)) +
  scale_x_log10() +
  ggtitle("Mod log.rich")

Residuals

resid <- residuals(grow2, type = "standardized")
resid <-  data.frame(resid[,1:4])
resid$log.rich <- diridatag$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")

ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
  facet_grid(~term)

MORT

predirim <-  mcomb %>% 
  filter(vital == "mort", data == "exclude") %>%
  select(fplot, term, VPC, richness.rarefaction,lat) 
diridatam <- predirim %>%
  pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
  mutate(log.rich = log(richness.rarefaction),
         log.rich.o=log.rich) %>%
  mutate_at(vars( log.rich), scale)

vpcm <- DR_data(diridatam[,c("quadrat", "sp", "quadrat:sp",  "Residual")])
#plot(vpc)
mort2 <- DirichReg(vpcm~log.rich,data=diridatam, model="alternative", base=4)
summary(mort2)
## Call:
## DirichReg(formula = vpcm ~ log.rich, data = diridatam, model = "alternative",
## base = 4)
## 
## Standardized Residuals:
##                 Min       1Q   Median      3Q     Max
## quadrat     -1.2846  -0.6061  -0.0454  0.4078  1.8812
## sp          -1.9682  -0.8456   0.1570  0.9680  2.6456
## quadrat:sp  -1.8027  -0.6593  -0.0102  0.2449  1.6414
## Residual    -2.9957  -0.2679   0.0598  0.6229  1.6857
## 
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.6901     0.1375 -12.295  < 2e-16 ***
## log.rich      0.4103     0.1473   2.786  0.00534 ** 
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.87149    0.10008  -8.708   <2e-16 ***
## log.rich    -0.14106    0.09865  -1.430    0.153    
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -1.2178     0.1135  -10.73   <2e-16 ***
## log.rich     -0.2471     0.1217   -2.03   0.0424 *  
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
## 
## PRECISION MODEL:
## ------------------------------------------------------------------
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.3788     0.1756   19.24   <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-likelihood: 88.64 on 7 df (43 BFGS + 1 NR Iterations)
## AIC: -163.3, BIC: -156
## Number of Observations: 21
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatam$log.rich),
                                      max(diridatam$log.rich), length.out=10))
pred <- predict(mort2, newdata = newdata, se=T)
confint(mort2)
## 
## 95% Confidence Intervals (original form)
## 
## - Beta-Parameters:
## Variable: quadrat
##                2.5%   Est.   97.5%
## (Intercept)  -1.960  -1.69  -1.421
## log.rich      0.122   0.41   0.699
## 
## Variable: sp
##                2.5%    Est.   97.5%
## (Intercept)  -1.068  -0.871  -0.675
## log.rich     -0.334  -0.141   0.052
## 
## Variable: quadrat:sp
##                2.5%    Est.   97.5%
## (Intercept)  -1.440  -1.218  -0.995
## log.rich     -0.486  -0.247  -0.009
## 
## Variable: Residual
##   variable omitted
## 
## - Gamma-Parameters
##              2.5%  Est.  97.5%
## (Intercept)  3.04  3.38   3.72
colnames(pred) <- colnames(vpcm)
newdata$log.rich.o <- newdata$log.rich*sd(diridatam$log.rich.o) +
  mean(diridatam$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredm <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")

newpredm %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
                                       "Residual")) %>%
  ggplot(aes(x=rich.o, y=pred, col=term)) +
  geom_line() + facet_grid(~term) +
  geom_point(data=predirim, aes(x=richness.rarefaction, y=VPC, col=term)) +
  scale_x_log10() +
  ggtitle("Mod log.rich")

Residuals

resid <- residuals(mort2, type = "standardized")
resid <-  data.frame(resid[,1:4])
resid$log.rich <- diridatam$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")

ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
  facet_grid(~term)

REC

predirir <-  mcomb %>% 
  filter(vital == "rec", data == "exclude") %>%
  select(fplot, term, VPC, richness.rarefaction,lat) 
diridatar <- predirir %>%
  pivot_wider(names_from = term, values_from = VPC) %>% ungroup() %>%
  mutate(log.rich = log(richness.rarefaction),
         log.rich.o=log.rich) %>%
  mutate_at(vars( log.rich), scale)

vpcr <- DR_data(diridatar[,c("quadrat", "sp", "quadrat:sp",  "Residual")])
#plot(vpc)
rec2 <- DirichReg(vpcr~log.rich,data=diridatar, model="alternative", base=4)
summary(rec2)
## Call:
## DirichReg(formula = vpcr ~ log.rich, data = diridatar, model = "alternative",
## base = 4)
## 
## Standardized Residuals:
##                 Min       1Q   Median      3Q     Max
## quadrat     -1.5124  -0.6320  -0.2358  0.8198  2.1508
## sp          -2.8666  -0.2882   0.1438  0.7867  1.6337
## quadrat:sp  -1.4254  -0.6122  -0.2100  0.3186  2.8831
## Residual    -1.5080  -0.6407  -0.0523  0.5100  3.1244
## 
## MEAN MODELS:
## ------------------------------------------------------------------
## Coefficients for variable no. 1: quadrat
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.586194   0.129600  -4.523 6.09e-06 ***
## log.rich     0.005999   0.140757   0.043    0.966    
## ------------------------------------------------------------------
## Coefficients for variable no. 2: sp
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.08768    0.11441  -0.766    0.443    
## log.rich    -0.71600    0.11995  -5.969 2.38e-09 ***
## ------------------------------------------------------------------
## Coefficients for variable no. 3: quadrat:sp
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.70344    0.13331  -5.277 1.32e-07 ***
## log.rich    -0.07571    0.14027  -0.540    0.589    
## ------------------------------------------------------------------
## Coefficients for variable no. 4: Residual
## - variable omitted (reference category) -
## ------------------------------------------------------------------
## 
## PRECISION MODEL:
## ------------------------------------------------------------------
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   3.2418     0.1825   17.77   <2e-16 ***
## ------------------------------------------------------------------
## Significance codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-likelihood: 71.27 on 7 df (43 BFGS + 1 NR Iterations)
## AIC: -128.5, BIC: -121.9
## Number of Observations: 19
## Links: Logit (Means) and Log (Precision)
## Parametrization: alternative
newdata <- expand.grid(log.rich = seq(min(diridatar$log.rich),
                                      max(diridatar$log.rich), length.out=10))
pred <- predict(rec2, newdata = newdata, se=T)
colnames(pred) <- colnames(vpcr)
newdata$log.rich.o <- newdata$log.rich*sd(diridatar$log.rich.o) +
  mean(diridatar$log.rich.o)
newdata$rich.o <- exp(newdata$log.rich.o)
newpredr <- cbind(newdata,pred) %>% pivot_longer(`quadrat`:`Residual`, names_to="term", values_to="pred")

newpredr %>% mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp",
                                       "Residual")) %>%
  ggplot(aes(x=rich.o, y=pred, col=term)) +
  geom_line() + facet_grid(~term) +
  geom_point(data=predirir, aes(x=richness.rarefaction, y=VPC, col=term)) +
  scale_x_log10() +
  ggtitle("Mod log.rich")

Residuals

resid <- residuals(rec2, type = "standardized")
resid <-  data.frame(resid[,1:4])
resid$log.rich <- diridatar$log.rich
resid <- resid %>% pivot_longer(1:4, names_to="term", values_to="resid")

ggplot(resid, aes(x=log.rich, y=resid)) + geom_point() +
  facet_grid(~term)

FIGURE

Calculated prediction interval

source(here("scripts/prediction_intervals_dirichlet_exclude.R"), local=T)
#load("results/prediction_intervals_dirichlet_exclude.Rdata") #  quants
predis <- bind_rows(list(grow=newpredg, mort=newpredm, rec=newpredr),
                    .id="vital") %>% 
  mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp","Residual"))

pontos <- bind_rows(list(grow=predirig, mort= predirim, rec=predirir), .id="vital")%>% 
  mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))
quants <- bind_rows(list(grow=quantg,mort=quantm, rec=quantr), .id="vital")%>% 
  mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))
# pvalues
res <- data.frame(vital = c("grow", "mort", "rec"),
                  P = NA,
                  term = "Residual")

test <- summary(grow2)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsg <- vals[ grep("log.rich", vals$coef),] %>% 
  mutate(term=c("quadrat", "sp", "quadrat:sp"))
test <- summary(mort2)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsm <- vals[ grep("log.rich", vals$coef),] %>% 
  mutate(term=c("quadrat", "sp", "quadrat:sp"))
test <- summary(rec2)
vals <- as.data.frame(test$coef.mat)[1:6,]
vals$coef <- names(test$coefficients)[1:6]
valsr <- vals[ grep("log.rich", vals$coef),] %>% 
  mutate(term=c("quadrat", "sp", "quadrat:sp")) 

pvals <- bind_rows(list(grow=valsg,mort=valsm,rec=valsr), .id="vital") %>%
  dplyr::select(vital,`Pr(>|z|)`, term) %>% rename(P = `Pr(>|z|)`)
pvals <- rbind(pvals,res) %>% 
  mutate(term = fct_relevel(term, "quadrat", "quadrat:sp", "sp", "Residual"))

pvals$x <- 300
pvals$y <- 0.70
pvals$sig <- paste0("p = ", round(pvals$P,3))
#pvals$sig[pvals$P >0.0166] <- "ns"
pvals$sig[pvals$term=="Residual"] <- ""
pvals$sig[pvals$sig=="p = 0"] <- "p < 0.0001"
library(wesanderson) #better colour palette
# Gradient color para latitutde
pal <- wes_palette("Zissou1",20, type = "continuous")[20:1] # azul frio-temperado
fdiri_lat <-ggplot(predis, aes(x=rich.o, y=pred)) +
  geom_line(size=1) +
  facet_grid(vital~term, labeller=grpvit) +
  geom_smooth(data=quants,aes(x=rich.o, y=lower), se=F, size=0.1)+
  geom_smooth(data=quants,aes(x=rich.o, y=upper), se=F, size=0.1)  +
  geom_ribbon(data=quants,aes(x=rich.o, ymin=lower, ymax=upper,
                              y=mean), alpha=0.05, fill="blue",
              size=0)+
  geom_point(data=pontos, aes(x=richness.rarefaction, y=VPC, fill=abs(lat)),
             col="black", size=2.5, pch=21)+
  scale_fill_gradientn(colors=pal) +
  scale_y_continuous(minor_breaks = NULL) +
  scale_x_log10() +
  geom_text(data=pvals, aes(x=x,y=y, label=sig), size=3, hjust=1,vjust=0,
            col="black")+
  
  xlab("Species richness (log10 scale)")+
  ylab("VPC") +
  theme_cowplot() +
  theme(panel.background = element_rect(colour="lightgray"),
        legend.position = "none")
fdiri_lat

png(here("figures/FIG_S5.5_dirichlet_regressions_lat_exclude.png"), height = 600, width=800, res=100)
fdiri_lat
dev.off()
## quartz_off_screen 
##                 2